
source("1.2.1.1_MLE_Point.R")
source("1.2.1.2_IF.R")
source("getProbScalarRD.R")

RDMLEst = function(param, y, x, va, vb, max.step = 1000, thres = 1e-6,
                   alpha.start = NULL, beta.start = NULL) {
    
    pa = dim(va)[2]; pb = dim(vb)[2]
    
    weights = rep(1, length(x))
    
    ## starting values for parameter optimization
    if (is.null(alpha.start)) 
        alpha.start = rep(0, pa)
    if (is.null(beta.start)) 
        beta.start = rep(0, pb)
    
    ### point estimate
    mle = max.likelihood(param, y, x, va, vb, alpha.start, beta.start, 
        weights, max.step, thres, pa, pb)
    point.est = mle$par
    converged = mle$convergence
    # print(point.est)
    alpha.ml = point.est[1:pa]
    beta.ml = point.est[(pa + 1):(pa + pb)]
    
    ### Computing Fisher Information:

    IF.alpha.beta = IF.alpha.beta(y, x, va, vb, vc, alpha.ml, beta.ml)
    
    return(list(point.est = c(alpha.ml,beta.ml), IF.alpha.beta = IF.alpha.beta))
    
} 
